Table of figures

Supporting	Main text
Figure S 1	Figure 1
Figure S 2	Figure 2
Figure S 3	Figure 3
Figure S 4	Figure 4
Figure S 5
Figure S 6
Figure S 7
Figure S 8
Figure S 9
Figure S 10
Figure S 11
Figure S 12

Table of tables

Unfortunately links only redirect to tabsets open by default but not hidden ones. Please look for tables related to analyses by phylum or by compartment in their respective tabsets.

Supporting Tables
Table S 1
Table S 2
Table S 3
Table S 4
Table S 5
Table S 6
Table S 7
Table S 8
Table S 9
Table S 10
Table S 11
Table S 12
Table S 13
Table S 14

Packages

Code

library(dplyr)
library(tidyr)
library(vegan)
library(ggplot2)
library(targets)
library(ggrepel)
library(patchwork)
theme_set(theme_light())

Macrofauna data availabity

Code

tar_load(faunadata)
knitr::kable(head(faunadata))

First lines of the original macrofaunal data
Point	Site	Replicate	Method	Habitat	Year	Season	Date	Species	Abundance
Belle-ile 1	Belle-ile	1	Grab	Maerl	2007	Spring	2007-03-11 00:00:00	Abludomelita gladiosa	1
Belle-ile 1	Belle-ile	1	Grab	Maerl	2007	Spring	2007-03-11 00:00:00	Abra alba	1
Belle-ile 1	Belle-ile	1	Grab	Maerl	2007	Spring	2007-03-11 00:00:00	Ampharete spp.	1
Belle-ile 1	Belle-ile	1	Grab	Maerl	2007	Spring	2007-03-11 00:00:00	Amphipholis squamata	2
Belle-ile 1	Belle-ile	1	Grab	Maerl	2007	Spring	2007-03-11 00:00:00	Anapagurus hyndmanni	1
Belle-ile 1	Belle-ile	1	Grab	Maerl	2007	Spring	2007-03-11 00:00:00	Aonides oxycephala	3

Code

samplespoint <- faunadata %>% 
  mutate(Date=as.factor(format.Date(faunadata$Date,"%Y"))) %>% 
  distinct(Site, Year, Point, Replicate) %>%
  group_by(Site, Year, Point) %>%
  mutate(Year = as.factor(Year)) %>% 
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(State = factor(case_when(count==3 ~ "Good", count < 3 ~ "Sample(s) missing", count > 3 ~ "Too many samples")))

ggplot(samplespoint, aes(x=Year,y=Point,fill=State)) +
  geom_tile(alpha=0.5) +
  scale_fill_manual(values=c("palegreen3","lightblue", "indianred")) +
  geom_text(aes(label=count)) +
  theme(text=element_text(size =15), axis.text.x = element_text(angle = 90, size = 15), legend.position = "bottom")

Supporting Figure 1: Number of grab samples taken every year at each point. Empty cells show points where no samples were collected at a given date due to field constraints.

Habitat Complexity

Take a look at the structure of the final complexity dataset

Code

tar_load(comp)
knitr::kable(head(comp))

First lines of the complexity dataset
Mini	Site	Point	Sample	Branching_density	Broken_density	D_bin	D_gray	L	I	S	Total_Density	Dry_weight	Sphericity	DR1	DR2	DR3
MBE_1	Belle-ile	Belle-ile 1	1	5.8	0.0	1.1231	2.6562	12.85333	9.143333	8.303333	181696	183.61	0.8371338	0.6460062	0.8153846	0.7113589
MBE_1	Belle-ile	Belle-ile 1	2	1.0	0.0	1.0654	2.6463	16.07667	7.883333	5.393333	181696	183.61	0.6122605	0.3354758	0.7669267	0.4903587
MBE_1	Belle-ile	Belle-ile 1	3	3.0	0.2	1.0722	2.6450	15.30000	8.103333	7.190000	181696	183.61	0.7470808	0.4699346	0.8873818	0.5296296
MBE_1	Belle-ile	Belle-ile 1	4	4.8	0.4	1.1078	2.6492	12.91667	8.476667	6.203333	181696	183.61	0.7057078	0.4802581	0.6613704	0.6562581
MBE_1	Belle-ile	Belle-ile 1	5	4.2	0.0	1.1305	2.7019	24.54667	14.523333	8.970000	181696	183.61	0.6088477	0.3654264	0.6434838	0.5916621
MBE_1	Belle-ile	Belle-ile 1	6	3.0	0.4	1.1198	2.6641	16.35667	9.416667	6.783333	181696	183.61	0.6684949	0.4147137	0.7249304	0.5757082

Most variables are not normally distributed so we apply a box-cox (Box and Cox 1964) transformation and check correlations between the transformed variables

Code

tar_read(pairscomp)

Correlogram of box-cox transformed HC variables

PCA of all observations

Let’s check it without removing highly correlated variables

Code

tar_read(pcatot)

PCA of all observations and all variables

Let’s check the biplots of the selected variables

Code

tar_read(pcasel)

Supporting Figure 2: PCA of all observations after variable selection

PCA of median values

Let’s check the biplots of the selected variables

Figure 1

Code

tar_read(pcamed)

Correlation between the centroids of the total PCA and the scores of the reduced PCA

Code

tar_load(corcomp)
tar_load(rv)
plot(corcomp[[1]])
plot(corcomp[[2]])
tar_read(cortestcomp1)
tar_read(cortestcomp2)

RV coefficient = 0.8744522, p-value = 6.8861509^{-12}

Cluster of median site values

Let’s check if we find similar groupings using a ward’s cluster

Code

tar_read(clustcomp)

Supporting Figure 4: Ward’s cluster of median complexity values at the site level

Check correlations between HC and other environmental variables

Code

tar_read(correl)

Supporting Figure 5: Correlogram of HC and physical environmental variables

Local macrofaunal diversity

Code

# Total richness
tar_load(faunaclass)
total_rich <- nrow(faunaclass %>% 
  distinct(., Species)) #725 species

# Total richness by sediment position

epi_rich <- nrow(faunaclass %>% 
  filter(Position == "Epifauna") %>% 
  distinct(., Species)) #341 epifaunal species


inf_rich <- nrow(faunaclass %>% 
  filter(Position == "Infauna") %>% 
  distinct(., Species)) #266  infaunal species


int_rich <- nrow(faunaclass %>% 
  filter(Position == "Interstice") %>% 
  distinct(., Species)) #118 interstitial species

There are 725 species in total, 341 being epifaunal, 266 being infaunal, and 118 interstitial.

Taxa from 14 phyla were identified, and Annelida, Arthropoda, Mollusca and Echinodermata were dominant, although Chordata, Cnidaria and Nemertea were also abundant at some sites. Some beds were consistently dominated by polychaetes, while others were always dominated by crustaceans, and a few were more even over the years with dominance shifts between Annelida, Arthropoda and Mollusca (Figure S 6). Maerl annelid communities were dominated by the mobile predators eunicids (Eunicidae, Lumbrineridae, Dorvilleidae), the crawling surface deposit feeders spionids (Spionidae), the tube dwelling suspension feeding sabellids (Serpullidae), and the subsurface deposit feeding knot worms (Polygordidae; Figure S 7 - A). As for arthropods (Figure S 7 - B), the decapods were dominant in most beds, namely the filter-feeding porcelain crabs Pisidia longicornis (Linnaeus, 1767; Porcellanidae), which was consistently dominant. The subsurface deposit feeding Corophiidae amphipods and Apseudiidae tanaids, as well as the scavenging Lysianassoidea amphipods were also abundant in most sites. Finally, maerl mollusc communities (Figure S 7 - C) were mostly dominated by small grazing and deposit feeding gastropods (Cerithiidae, Rissoidea and Trochidae), by endo and epifaunal filter-feeding pectinid bivalves (namely Veneridae), and by chitons (Leptochitonidae).

Code

tar_read(relab_phy)
tar_read(relab_ord)

(p.s.: it is normal to have NAs for some orders as higher classification for some of the species is ambiguous or unknown.)

Code

tar_read(relab_arth)
tar_read(relab_ann) 
tar_read(relab_mol)

Effects of HC in local diversity

Figure 2

Code

tar_read(figure3) #figure 2-a
tar_read(psdens) #figure 2-b

Figure 2 (main text)

Observed richness

Check if the relationships between HC and species richness are the same for the total macrofauna as well as for the main phyla and the main compartments (different positions in sediment) separately. Main results shown in Figure 2

Now we first forward select environmental variables and then add sites as a random factor to control for pseudo-replication and other site-related variables that were not taken into account.

Code

tar_read(totrichpc1)

OLS linear regression of Total Macrofauna Richness as a function of Rhodolith complexity (Bed complexity is shown in Figure 2)

We check models with quadratic functions of both PC1 and PC2 to test for a unimodal relationship. Unimodal trends are never significant, suggesting an overall linear relationship. A significant interaction is found between PC1² and both PC2 and PC2² in the last model in Table S 1. Since there’s a very slight difference in AICs and BICs between the last two models, we keep the more parsimonious one.

Code

tar_load(c(modfixrich1, modfixrichpoly1, modfixrich2, modfixrichpoly2, modfixrich12, modfixrichpoly12))

jtools::export_summs(
  modfixrich1, 
  modfixrichpoly1,
  modfixrich2,
  modfixrichpoly2,
  modfixrich12,
  modfixrichpoly12,
  model.names = c("S ~ PC1", "S ~ PC1+PC1^2", "S ~ PC2", "S ~ PC2+PC2^2", "S ~ PC1*PC2","S ~ (PC1+PC1^2)*(PC2+PC2^2)"),
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "Rhodolith \n complexity \n (PC1)" = "PC1_score",
    "Rhodolith \n complexity \n (PC1)" = "`poly(PC1_score, 2)`1",
    "PC1^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "Bed complexity (PC2)" = "`poly(PC2_score, 2)`1",
    "PC2^2" = "`poly(PC2_score, 2)`2",
    "PC1:PC2" = "PC1_score:PC2_score",
    "PC1^2:PC2" ="`poly(PC1_score, 2)`2:`poly(PC2_score, 2)`1",
    "PC1:PC2^2" ="`poly(PC1_score, 2)`1:`poly(PC2_score, 2)`2",
    "PC1^2:PC2^2" ="`poly(PC1_score, 2)`2:`poly(PC2_score, 2)`2"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 1: Model statistics of the linear and unimodal models of Richness
	S ~ PC1		S ~ PC1+PC1^2		S ~ PC2		S ~ PC2+PC2^2		S ~ PC1*PC2		S ~ (PC1+PC1^2)*(PC2+PC2^2)
Rhodolith complexity (PC1)	8.22 ***	(1.00)	8.21 ***	(1.00)					8.38 ***	(0.91)	6.96 ***	(1.02)
PC1^2			1.37	(1.00)							-1.25	(1.06)
Bed complexity (PC2)					8.14 ***	(1.00)	8.14 ***	(1.00)	8.03 ***	(0.92)	8.89 ***	(1.17)
PC2^2							0.44	(1.00)			1.30	(1.20)
PC1:PC2									-0.87	(0.97)
PC1^2:PC2											4.02 **	(1.25)
PC1:PC2^2											-1.80	(1.30)
PC1^2:PC2^2											-3.58 **	(1.19)
N	348		348		348		348		348		348
DF	346.00		345.00		346.00		345.00		344.00		339.00
AIC	3028.54		3028.65		3029.73		3031.53		2956.83		2945.82
BIC	3040.09		3044.06		3041.28		3046.94		2976.10		2984.35
Adjusted R2	0.16		0.16		0.16		0.16		0.32		0.35
F	67.23		34.64		65.82		32.93		55.56		24.45
p-value	0.00		0.00		0.00		0.00		0.00		0.00
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_load(c(modfixrich12, modfrsel, modrrsel))

jtools::export_summs(
  modfixrich12,
  modfrsel,
  modrrsel,
  model.names = c("S ~ Complexity", "S ~ Complexity + Environment", "S ~ Complexity + Environment + (1|Site)"),
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Depth",
    "Exposure" = "Fetch_max",
    "Mean bottom temperature" = "T_mean",
    "Organic Matter %" = "OM",
    "Year"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 2: Model statistics of the final selected models of Richness
	S ~ Complexity		S ~ Complexity + Environment		S ~ Complexity + Environment + (1\|Site)
Rhodolith complexity (PC1)	8.38 ***	(0.91)	2.89 *	(1.20)	2.08	(1.30)
Bed complexity (PC2)	8.03 ***	(0.92)	4.30 ***	(1.15)	4.05 *	(1.58)
PC1:PC2	-0.87	(0.97)	0.04	(0.91)	2.04	(1.10)
Depth			-5.33 ***	(1.25)	-3.06	(2.45)
Exposure			-2.75	(1.41)	-1.54	(2.83)
Mean bottom temperature			2.05 *	(0.88)	1.38	(0.86)
Organic Matter %			2.56 *	(1.12)	3.72 **	(1.23)
Year			1.71	(0.87)	1.78 *	(0.83)
N	348		348		348
DF	344.00		339.00		337.00
AIC	2956.83		2905.35		2872.96
BIC	2976.10		2943.88		2915.33
Marginal R2					0.29
Conditional R2	0.33		0.44		0.43
Adjusted R2	0.32		0.42
F	55.56		32.69
p-value	0.00		0.00
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_load(c(modrrann, modrrart, modrrmol, modrrepi, modrrinf, modrrint))

jtools::export_summs(
  modrrann,
  modrrart,
  modrrmol,
  model.names = c("Annelida", "Arthropoda", "Mollusca"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 3: Model statistics of richness as a function of HC for each main phyla
	Annelida		Arthropoda		Mollusca
Rhodolith complexity (PC1)	1.50 *	(0.66)	0.43	(0.44)	-0.15	(0.48)
Bed complexity (PC2)	2.50 **	(0.82)	0.90	(0.54)	0.11	(0.60)
PC1:PC2	1.50 **	(0.56)	0.59	(0.36)	0.06	(0.40)
N	348		348		348
Residual DF	342.00		342.00		342.00
Marginal R2	0.11		0.02		0.00
Conditional R2	0.38		0.52		0.51
p-value
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(annrichpc1)
tar_read(annrichpc2) 
tar_read(artrichpc1)
tar_read(artrichpc2)
tar_read(molrichpc1)
tar_read(molrichpc2)

Main phyla richness as a function of HC

Code

jtools::export_summs(
  modrrepi,
  modrrinf,
  modrrint,
  model.names = c("Epifauna", "Infauna", "Interstitial fauna"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
    coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 4: Model statistics of richness as a function of HC for each compartiment
	Epifauna		Infauna		Interstitial fauna
Rhodolith complexity (PC1)	0.78	(0.74)	0.37	(0.56)	0.80 **	(0.29)
Bed complexity (PC2)	1.47	(0.93)	1.35	(0.70)	1.35 ***	(0.35)
PC1:PC2	1.40 *	(0.62)	0.76	(0.47)	0.36	(0.25)
N	348		348		347
Residual DF	342.00		342.00		341.00
Marginal R2	0.03		0.03		0.15
Conditional R2	0.52		0.44		0.31
p-value
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(epirichpc1)
tar_read(epirichpc2) 
tar_read(infrichpc1)
tar_read(infrichpc2)
tar_read(intrichpc1)
tar_read(intrichpc2)

Species richness as a function of HC at each faunal compartment

As we can see, the marginal R2 for the models containing only a subset of the macrofauna are quite low. However, although the effects of HC are less evident in these cases, The pattern for most subsets is very similar to the general pattern, with very slight positive effect of rhodolith complexity and a more important positive effect of bed complexity. Nevertheless Rhodolith complexity seems to play a more important role in driving annelid richness than the other main phyla.

Density of macrofauna

Code

tar_load(c(modfixdens12, modfdsel, modrdsel))
jtools::export_summs(
  modfixdens12,
  modfdsel,
  modrdsel,
  model.names = c("Density ~ Complexity", "Density ~ Complexity + Environment", "Density ~ Complexity + Environment + (1|Site)"),
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
    coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Mean current velocity" = "Current_mean",
    "Depth",
    "Exposure" = "Fetch_max",
    "Gravel",
    "Organic Matter %" = "OM",
    "Mean Temperature" = "T_mean",
    "Year"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 5: Model statistics of macrofaunal density as a function of HC
	Density ~ Complexity		Density ~ Complexity + Environment		Density ~ Complexity + Environment + (1\|Site)
Rhodolith complexity (PC1)	0.42 ***	(0.05)	0.21 ***	(0.06)	0.06	(0.05)
Bed complexity (PC2)	0.09	(0.05)	-0.07	(0.05)	0.13	(0.07)
PC1:PC2	-0.31 ***	(0.05)	-0.14 **	(0.05)	0.04	(0.04)
Mean current velocity			-0.27 ***	(0.05)	-0.50	(0.24)
Depth			-0.29 ***	(0.06)	0.08	(0.14)
Exposure			0.20 **	(0.07)	0.18	(0.18)
Gravel			0.24 ***	(0.04)	0.05	(0.05)
Organic Matter %			0.20 ***	(0.05)	0.14 **	(0.05)
Mean Temperature			0.11 **	(0.04)	0.06	(0.03)
Year			0.03	(0.04)	0.05	(0.03)
N	348		348		348
DF	344.00		337.00		335.00
AIC	889.58		767.68		687.22
BIC	908.84		813.91		737.30
Marginal R2					0.24
Conditional R2	0.26		0.50		0.69
Adjusted R2	0.25		0.48
F	39.42		33.22
p-value	0.00		0.00
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_load(c(modfdann, modrdann, modfdart, modrdart, modfdmol, modrdmol, modfdepi, modrdepi, modfdinf, modrdinf, modfdint, modrdint))

jtools::export_summs(
  modfdann,
  modrdann,
  modfdart,
  modrdart,
  modfdmol,
  modrdmol,
  model.names = c("Annelida LM","Annelida LMM", "Arthropoda LM", "Arthropoda LMM", "Mollusca LM", "Mollusca LMM"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Mean current velocity" = "Current_mean",
    "Depth",
    "Exposure" = "Fetch_max",
    "Gravel",
    "Organic Matter %" = "OM",
    "Mean Temperature" = "T_mean",
    "Year"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 6: Model statistics of Density as a function of HC for each main phyla
	Annelida LM		Annelida LMM		Arthropoda LM		Arthropoda LMM		Mollusca LM		Mollusca LMM
Rhodolith complexity (PC1)	0.05	(0.05)	0.07	(0.05)	0.31 ***	(0.07)	0.03	(0.05)	0.02	(0.05)	-0.02	(0.06)
Bed complexity (PC2)	0.11 *	(0.05)	0.11	(0.07)	0.00	(0.06)	0.10	(0.06)	0.31 ***	(0.05)	0.18 **	(0.07)
PC1:PC2	-0.01	(0.04)	0.05	(0.05)	-0.18 ***	(0.05)	0.06	(0.04)	-0.05	(0.04)	-0.01	(0.05)
Mean current velocity					-0.29 ***	(0.06)	-0.10	(0.22)
Depth	-0.31 ***	(0.05)	-0.21 *	(0.10)					-0.42 ***	(0.06)	-0.36 ***	(0.10)
Exposure	-0.29 ***	(0.06)	-0.37 **	(0.12)	0.30 ***	(0.07)	0.31	(0.18)
Gravel					0.31 ***	(0.05)	0.02	(0.04)	0.17 ***	(0.04)	0.05	(0.05)
Organic Matter %	0.06	(0.05)			0.22 ***	(0.06)	0.10 *	(0.05)	0.09 *	(0.05)	0.12 *	(0.06)
Mean Temperature	0.09 *	(0.04)	0.07 *	(0.04)					0.10 **	(0.04)	0.09 *	(0.04)
Year	0.15 ***	(0.04)	0.16 ***	(0.03)	0.06	(0.04)	0.07 *	(0.03)	-0.06	(0.04)	-0.06	(0.04)
N	348		348		348		348		348		348
Residual DF	338.00		338.00		339.00		337.00		339.00		337.00
Marginal R2			0.47				0.04				0.35
Conditional R2	0.59		0.58		0.37		0.75		0.54		0.47
Adjusted R2	0.58				0.36				0.53
F	53.91				24.94				49.88
p-value	0.00				0.00				0.00
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(anndenspc1)
tar_read(anndenspc2) 
tar_read(artdenspc1)
tar_read(artdenspc2)
tar_read(moldenspc1)
tar_read(moldenspc2)

Main phyla density as a function of HC

Code

jtools::export_summs(
  modfdepi,
  modrdepi,
  modfdinf,
  modrdinf,
  modfdint,
  modrdint,
  model.names = c("Epifauna LM", "Epifauna LMM", "Infauna LM", "Infauna LMM", "Interstitial fauna LM", "Interstitial fauna LMM"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Mean current velocity" = "Current_mean",
    "Depth",
    "Exposure" = "Fetch_max",
    "Gravel",
    "Organic Matter %" = "OM",
    "Mud",
    "Mean Temperature" = "T_mean",
    "Temperature Variation (sd)" = "T_sd",
    "Year"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 7: Model statistics of density of macrofauna as a function of HC for each compartiment
	Epifauna LM		Epifauna LMM		Infauna LM		Infauna LMM		Interstitial fauna LM		Interstitial fauna LMM
Rhodolith complexity (PC1)	0.31 ***	(0.06)	0.07	(0.05)	0.10	(0.05)	0.01	(0.06)	0.06	(0.05)	0.04	(0.05)
Bed complexity (PC2)	0.01	(0.05)	0.11	(0.07)	0.17 ***	(0.05)	0.14	(0.07)	0.22 ***	(0.04)	0.20 **	(0.06)
PC1:PC2	-0.09	(0.05)	0.07	(0.04)	-0.11 *	(0.05)	-0.00	(0.05)	-0.09 *	(0.04)	-0.03	(0.05)
Mean current velocity	-0.28 ***	(0.06)	-0.23	(0.21)
Depth					-0.48 ***	(0.06)	-0.20	(0.12)
Exposure	0.27 ***	(0.07)	0.21	(0.18)					-0.34 ***	(0.06)	-0.40 ***	(0.09)
Gravel	0.38 ***	(0.05)	0.07	(0.05)	-0.25 ***	(0.05)	-0.10	(0.05)	0.26 ***	(0.04)	0.19 ***	(0.05)
Organic Matter %	0.27 ***	(0.05)	0.16 **	(0.05)					0.03	(0.05)	0.04	(0.05)
Mud									0.33 ***	(0.05)	0.25 ***	(0.06)
Mean Temperature	0.08	(0.04)	0.04	(0.03)	0.08	(0.04)	0.05	(0.04)
Temperature Variation (sd)					-0.01	(0.04)	-0.01	(0.06)
Year	0.00	(0.04)	0.02	(0.03)	0.09 *	(0.04)	0.09 *	(0.04)	0.14 ***	(0.03)	0.14 ***	(0.03)
N	348		348		348		348		348		348
Residual DF	338.00		336.00		339.00		337.00		339.00		337.00
Marginal R2			0.09				0.11				0.56
Conditional R2	0.42		0.66		0.43		0.49		0.60		0.62
F	27.06				32.32				64.41
p-value	0.00				0.00				0.00
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(epidenspc1)
tar_read(epidenspc2) 
tar_read(infdenspc1)
tar_read(infdenspc2)
tar_read(intdenspc1)
tar_read(intdenspc2)

Density of macrofauna as a function of HC at each compartment

Code

tar_read(psdtrait)

Macrofaunal densities as a function of HC for each faunal compartment

N2

Code

tar_read(psN2)

Supporting Figure 8: Regression coefficients of LMs and LMM of Simpson’s inverse diversity (N2) as a function of HC and physical environmental variables.

J

Code

tar_read(psJ)

Supporting Figure 9: Regression coefficients of LMs and LMM of Pielou’s evenness (J) as a function of HC and physical

Regional patterns of biodiversity

PCA

The first axis of the PCA of the log-chord transformed community data (Figure S 10) highlights a very similar gradient in macrofaunal diversity to that of rhodolith complexity (Figure 1) and accounts for around 14 % of the variance in assemblages. The second axis mostly captures the latitudinal gradient in the study area, which explains around 9.3 % of the total variance. These gradients are mostly explained by differences in the presence and abundance of a few key species (the ones with at least 55 % of their variance represented in the first two axes). The crustacean Athanas nitescens (Leach, 1814) is very rare or never found in most sites except for Rozegat, Keraliou and Trevignon. These sites are also grouped by the presence of Pisidia longicornis, which is found in most sites but is a lot more abundant at the bay of Brest and Trevignon, reaching up to 90 % of relative abundance during some years in the latter. The crustacean Leptocheirus tricristatus (Chevreux, 1887) is present in most sites, with the exception of Rozegat and Keraliou, and is more abundant in Belle-Ile and Meaban (reaching up to 12 % of relative abundance). Finally, the polychaete Glycera lapidum is more abundant in Saint-Brieuc and Camaret but is very rare or absent in sites at the other end of the gradient.

Code

tar_read(faunapca)

Supporting Figure 10: Principal Component Analysis of log-chord transformed macrofaunal species abundances. The first two axes represent around 23 % of the variance in community composition and structure. Left panel: distance biplot (scaling one) with outer plots representing the density of points for each site. Right panel: correlation biplot (scaling 2) showing species with for which at least 0.55 of their total variance is represented in the ordination plane.

Cluster

Code

tar_read(faunaclust)

Ward’s cluster of median complexity values at the site level

Effects of HC on regional diversity patterns

Redundancy analysis

After variable selection, we add sites as a factor to understand if we describe them well with the environmental variables we chose. We can look at the first two axes to have an idea of the relationships between the different explanatory variables.

Code

tar_read(rdafselcomp)

          variables order         R2      R2Cum  AdjR2Cum        F pvalue
1        Sphericity     4 0.07252058 0.07252058 0.0698400 27.05410  0.001
2 Branching_density     1 0.05583365 0.12835423 0.1233012 22.09912  0.001
3     Total_Density     3 0.03274661 0.16110084 0.1537849 13.42811  0.001
4                 L     2 0.02417722 0.18527806 0.1757769 10.17867  0.001
5               DR3     5 0.01151774 0.19679579 0.1850530  4.90419  0.001

Code

tar_read(rdafseltemp)

  variables order         R2      R2Cum   AdjR2Cum         F pvalue
1      T_sd     2 0.03456694 0.03456694 0.03177668 12.388392  0.001
2    T_mean     1 0.01675283 0.05131977 0.04582018  6.092385  0.001

Code

tar_read(rdafselgranulo)

  variables order         R2     R2Cum   AdjR2Cum        F pvalue
1       Mud     1 0.08700660 0.0870066 0.08436789 32.97317  0.001
2    Gravel     2 0.04160464 0.1286112 0.12355971 16.47210  0.001
3        OM     3 0.01247367 0.1410849 0.13359437  4.99577  0.001

Code

tar_read(rdafselhydro)

     variables order         R2      R2Cum   AdjR2Cum        F pvalue
1 Current_mean     1 0.06103081 0.06103081 0.05831702 22.48919  0.001
2    Fetch_max     2 0.04553573 0.10656654 0.10138721 17.58366  0.001

Code

tar_load(rdasite)
tar_read(triplotrda)

RDA triplot showing the variability in species composition related to physical environmental variables and habitat complexity. The first two axes are shown.

Adjusted R² `r RsquareAdj(rdasite)$adj.r.squared We check the model’s validity. Only 99 permutations were chosen for the validation by axis and by term here as they are really time consuming but the results are the same with 999 permutations. (Output is collapsed, click on code to see)

Code

tar_read(aovsite)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + Total_Density + L + OM + Year + Site, data = envcomp)
##           Df Variance      F Pr(>F)    
## Model     23  0.28507 12.289  0.001 ***
## Residual 324  0.32679                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tar_read(aovaxis)
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 99
## 
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + Total_Density + L + OM + Year + Site, data = envcomp)
##           Df Variance       F Pr(>F)   
## RDA1       1  0.07620 75.5463   0.01 **
## RDA2       1  0.04629 45.8898   0.01 **
## RDA3       1  0.04244 42.0730   0.01 **
## RDA4       1  0.02841 28.1713   0.01 **
## RDA5       1  0.01876 18.5971   0.01 **
## RDA6       1  0.01339 13.2736   0.01 **
## RDA7       1  0.01225 12.1454   0.01 **
## RDA8       1  0.00957  9.4930   0.01 **
## RDA9       1  0.00781  7.7480   0.01 **
## RDA10      1  0.00710  7.0381   0.01 **
## RDA11      1  0.00557  5.5256   0.01 **
## RDA12      1  0.00307  3.0402   0.01 **
## RDA13      1  0.00239  2.3719   0.04 * 
## RDA14      1  0.00202  2.0023   0.12   
## RDA15      1  0.00176  1.7401   0.38   
## RDA16      1  0.00154  1.5284   0.72   
## RDA17      1  0.00133  1.3138   0.96   
## RDA18      1  0.00113  1.1189   1.00   
## RDA19      1  0.00103  1.0169   1.00   
## RDA20      1  0.00088  0.8709   1.00   
## RDA21      1  0.00077  0.7609   1.00   
## RDA22      1  0.00073  0.7236   1.00   
## RDA23      1  0.00065  0.6482   1.00   
## Residual 324  0.32679                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tar_read(aovterm)
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 99
## 
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + Total_Density + L + OM + Year + Site, data = envcomp)
##                    Df Variance       F Pr(>F)   
## Mud                 1  0.05324 52.7815   0.01 **
## Depth               1  0.03189 31.6154   0.01 **
## Current_mean        1  0.03207 31.7934   0.01 **
## Branching_density   1  0.02296 22.7680   0.01 **
## Fetch_max           1  0.01678 16.6339   0.01 **
## T_sd                1  0.01202 11.9166   0.01 **
## Sphericity          1  0.00898  8.9014   0.01 **
## T_mean              1  0.00709  7.0263   0.01 **
## DR3                 1  0.00610  6.0457   0.01 **
## Gravel              1  0.00858  8.5018   0.01 **
## Total_Density       1  0.00236  2.3400   0.01 **
## L                   1  0.00491  4.8632   0.01 **
## OM                  1  0.00426  4.2272   0.01 **
## Year                1  0.00573  5.6773   0.01 **
## Site                9  0.06813  7.5051   0.01 **
## Residual          324  0.32679                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hierarchical partitioning

Figure 3

Code

tar_read(upsetmain)

Code

tar_read(upsetsep)

Supporting Figure 11: Hierarchical partitioning separating Rhodolith and Bed complexity

Community stability

Figure 4

Code

tar_read(bdplot)

Code

tar_read(bdtrait)
tar_read(bdphyl)

Effects of HC on community stability by compartiment

Effects of HC on community stability by phylum - not shown in main text

Code

tar_read(replplot)

Supporting Figure 12: Relative importance of the Replacement component of temporal BD

Effects of HC on community stability

Code

tar_load(c(modbdtot, modrdftot, modrptot))
jtools::export_summs(
  modbdtot,
  modrptot,
  modrdftot,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 8: Models statistics of temporal diversity metrics as a function of HC
	BDtotal		Replacement		Richness Difference
Rhodolith complexity (PC1)	-0.51 **	(0.16)	-0.42 *	(0.17)	-0.29	(0.18)
Rhodolith complexity (PC1)^2	-0.23	(0.16)	-0.10	(0.17)	-0.27	(0.19)
Bed complexity (PC2)	-0.18	(0.17)	-0.21	(0.19)	-0.02	(0.20)
PC1:PC2	0.28	(0.19)	0.11	(0.21)	0.32	(0.22)
PC1ˆ2:PC2	-0.07	(0.20)	-0.27	(0.22)	0.24	(0.24)
N	30		30		30
Residual DF	24.00		24.00		24.00
R2	0.43		0.32		0.20
Adjusted R2	0.31		0.18		0.03
F	3.55		2.23		1.20
p-value	0.02		0.08		0.34
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(bdtotpc1)
tar_read(bdtotpc2)
tar_read(rdtotpc1)
tar_read(rdtotpc2)
tar_read(rptotpc1)
tar_read(rptotpc2)

Code

tar_load(c(modbdann, modrdfann, modrpann))
jtools::export_summs(
  modbdann,
  modrpann,
  modrdfann,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 9: Models statistics of annelid temporal diversity metrics as a function of HC
	BDtotal		Replacement		Richness Difference
Rhodolith complexity (PC1)	-0.61 ***	(0.14)	-0.64 ***	(0.15)	-0.16	(0.19)
Rhodolith complexity (PC1)^2	-0.27	(0.15)	-0.11	(0.15)	-0.30	(0.19)
Bed complexity (PC2)	-0.08	(0.16)	0.01	(0.16)	-0.14	(0.20)
PC1:PC2	0.22	(0.17)	0.10	(0.18)	0.22	(0.23)
PC1ˆ2:PC2	-0.10	(0.19)	-0.22	(0.19)	0.12	(0.24)
N	30		30		30
Residual DF	24.00		24.00		24.00
R2	0.51		0.49		0.18
Adjusted R2	0.41		0.38		0.01
F	4.98		4.62		1.03
p-value	0.00		0.00		0.42
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(bdannpc1)
tar_read(bdannpc2)
tar_read(rdannpc1)
tar_read(rdannpc2)
tar_read(rpannpc1)
tar_read(rpannpc2)

Effects of HC on annelid community stability

Code

tar_load(c(modbdart, modrdfart, modrpart))
jtools::export_summs(
  modbdart,
  modrpart,
  modrdfart,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 10: Models statistics of arthropod temporal diversity metrics as a function of HC
	BDtotal		Replacement		Richness Difference
Rhodolith complexity (PC1)	-0.42 *	(0.17)	-0.29	(0.16)	-0.38 *	(0.18)
Rhodolith complexity (PC1)^2	-0.09	(0.18)	-0.20	(0.16)	0.03	(0.18)
Bed complexity (PC2)	-0.23	(0.19)	-0.40 *	(0.18)	-0.02	(0.20)
PC1:PC2	0.22	(0.21)	-0.13	(0.19)	0.41	(0.22)
PC1ˆ2:PC2	-0.11	(0.22)	-0.41	(0.21)	0.16	(0.23)
N	30		30		30
Residual DF	24.00		24.00		24.00
R2	0.31		0.39		0.25
Adjusted R2	0.16		0.27		0.09
F	2.14		3.11		1.57
p-value	0.10		0.03		0.21
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(bdartpc1)
tar_read(bdartpc2)
tar_read(rdartpc1)
tar_read(rdartpc2)
tar_read(rpartpc1)
tar_read(rpartpc2)

Effects of HC on arthropod community stability

Code

tar_load(c(modbdmol, modrdfmol, modrpmol))
jtools::export_summs(
  modbdmol,
  modrpmol,
  modrdfmol,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 11: Models statistics of mollusc temporal diversity metrics as a function of HC
	BDtotal		Replacement		Richness Difference
Rhodolith complexity (PC1)	-0.22	(0.16)	-0.06	(0.19)	-0.24	(0.18)
Rhodolith complexity (PC1)^2	-0.32	(0.16)	-0.16	(0.20)	-0.23	(0.19)
Bed complexity (PC2)	-0.35	(0.17)	-0.22	(0.21)	-0.20	(0.20)
PC1:PC2	0.30	(0.19)	0.19	(0.23)	0.15	(0.22)
PC1ˆ2:PC2	0.15	(0.21)	0.02	(0.25)	0.20	(0.24)
N	30		30		30
Residual DF	24.00		24.00		24.00
R2	0.40		0.13		0.19
Adjusted R2	0.27		-0.05		0.03
F	3.18		0.72		1.16
p-value	0.02		0.62		0.36
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(bdmolpc1)
tar_read(bdmolpc2)
tar_read(rdmolpc1)
tar_read(rdmolpc2)
tar_read(rpmolpc1)
tar_read(rpmolpc2)

Effects of HC on mollusc community stability

Code

tar_load(c(modbdepi, modrdfepi, modrpepi))
jtools::export_summs(
  modbdepi,
  modrpepi,
  modrdfepi,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 12: Models statistics of epifaunal temporal diversity metrics as a function of HC
	BDtotal		Replacement		Richness Difference
Rhodolith complexity (PC1)	-0.42 *	(0.18)	-0.25	(0.19)	-0.38 *	(0.18)
Rhodolith complexity (PC1)^2	-0.16	(0.18)	-0.06	(0.19)	-0.19	(0.19)
Bed complexity (PC2)	-0.04	(0.19)	-0.07	(0.21)	0.02	(0.20)
PC1:PC2	0.29	(0.21)	0.08	(0.23)	0.36	(0.22)
PC1ˆ2:PC2	-0.05	(0.23)	-0.26	(0.25)	0.18	(0.24)
N	30		30		30
Residual DF	24.00		24.00		24.00
R2	0.27		0.15		0.22
Adjusted R2	0.12		-0.03		0.06
F	1.82		0.82		1.38
p-value	0.15		0.55		0.27
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(bdepipc1)
tar_read(bdepipc2)
tar_read(rdepipc1)
tar_read(rdepipc2)
tar_read(rpepipc1)
tar_read(rpepipc2)

Effects of HC on epifaunal community stability

Code

tar_load(c(modbdinf, modrdfinf, modrpinf))
jtools::export_summs(
  modbdinf,
  modrpinf,
  modrdfinf,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 13: Models statistics of infaunal temporal diversity metrics as a function of HC
	BDtotal		Replacement		Richness Difference
Rhodolith complexity (PC1)	-0.56 ***	(0.14)	-0.43 *	(0.17)	-0.28	(0.18)
Rhodolith complexity (PC1)^2	-0.35 *	(0.14)	-0.25	(0.17)	-0.20	(0.19)
Bed complexity (PC2)	-0.26	(0.15)	-0.16	(0.18)	-0.19	(0.20)
PC1:PC2	0.18	(0.16)	0.09	(0.20)	0.17	(0.22)
PC1ˆ2:PC2	-0.09	(0.18)	-0.18	(0.22)	0.13	(0.24)
N	30		30		30
Residual DF	24.00		24.00		24.00
R2	0.57		0.33		0.19
Adjusted R2	0.48		0.20		0.03
F	6.30		2.41		1.15
p-value	0.00		0.07		0.36
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(bdinfpc1)
tar_read(bdinfpc2)
tar_read(rdinfpc1)
tar_read(rdinfpc2)
tar_read(rpinfpc1)
tar_read(rpinfpc2)

Effects of HC on infaunal community stability

Code

tar_load(c(modbdint, modrdfint, modrpint))
jtools::export_summs(
  modbdint,
  modrpint,
  modrdfint,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)

Supporting Table 14: Models statistics of interstitial fauna temporal diversity metrics as a function of HC
	BDtotal		Replacement		Richness Difference
Rhodolith complexity (PC1)	-0.33	(0.16)	-0.46 *	(0.17)	0.12	(0.18)
Rhodolith complexity (PC1)^2	-0.12	(0.17)	0.17	(0.17)	-0.31	(0.18)
Bed complexity (PC2)	-0.40 *	(0.18)	-0.26	(0.18)	-0.16	(0.20)
PC1:PC2	0.26	(0.20)	0.01	(0.20)	0.27	(0.22)
PC1ˆ2:PC2	-0.02	(0.21)	-0.32	(0.22)	0.30	(0.23)
N	30		30		30
Residual DF	24.00		24.00		24.00
R2	0.38		0.35		0.25
Adjusted R2	0.25		0.21		0.09
F	2.92		2.55		1.57
p-value	0.03		0.05		0.21
All continuous variables are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

tar_read(bdintpc1)
tar_read(bdintpc2)
tar_read(rdintpc1)
tar_read(rdintpc2)
tar_read(rpintpc1)
tar_read(rpintpc2)

Effects of HC on intersitial community stability

Jardim et al. (2023): Supporting Information

Table of figures

Table of tables

Packages

Macrofauna data availabity

Habitat Complexity

PCA of all observations

PCA of median values

Correlation between the centroids of the total PCA and the scores of the reduced PCA

Cluster of median site values

Check correlations between HC and other environmental variables

Local macrofaunal diversity

Effects of HC in local diversity

Observed richness

Density of macrofauna

N2

J

Regional patterns of biodiversity

PCA

Cluster

Effects of HC on regional diversity patterns

Redundancy analysis

Hierarchical partitioning

Community stability

Effects of HC on community stability

References